# Replication package for 
# Making Sense of Violence in Semi-Technologized Conventional Civil War
# Publication Tables II Replication Files
# Analyses run using RStudio Version: 2022.07.1+554

rm(list = ls())

######################################
### Load packages
######################################

library(psych)
library(mgcv)

# Set working directory to folder with replication files
setwd("Place path to data files here")

#################################################
# Use Boshin War data
#################################################

boshin_war <- read.csv("boshin_war.csv", na = " ")

# Table 1: Descriptive statistics
c <- describe(boshin_war, na.rm = FALSE, fast=TRUE)
print(c, digits = 3)

# Table 2: Zero-inflated Poisson estimations with thin-plate regression splines for coordinates
gam.model1 <- gam(vtype123 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=ziP(theta = NULL, link = "identity", b=0), data=boshin_war, method="GCV.Cp")
summary(gam.model1)
gam.model2 <- gam(vtype1 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=ziP(theta = NULL, link = "identity", b=0), data=boshin_war, method="GCV.Cp")
summary(gam.model2)
gam.model3 <- gam(vtype23 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=ziP(theta = NULL, link = "identity", b=0), data=boshin_war, method="GCV.Cp")
summary(gam.model3)
gam.model4 <- gam(vtype2 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=ziP(theta = NULL, link = "identity", b=0), data=boshin_war, method="GCV.Cp")
summary(gam.model4)
gam.model5 <- gam(vtype3 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=ziP(theta = NULL, link = "identity", b=0), data=boshin_war, method="GCV.Cp")
summary(gam.model5)
gam.model6 <- gam(vtype4 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + frontline + s(longitude, latitude), family=ziP(theta = NULL, link = "identity", b=0), data=boshin_war, method="GCV.Cp")
summary(gam.model6)

# Table 3: Logit estimations with thin-plate splines for coordinates
gam.model7 <- gam(dum_vtype123 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=binomial(link="logit"), data=boshin_war, method="GCV.Cp")
summary(gam.model7)
gam.model8 <- gam(dum_vtype1 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=binomial(link="logit"), data=boshin_war, method="GCV.Cp")
summary(gam.model8)
gam.model9 <- gam(dum_vtype23 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=binomial(link="logit"), data=boshin_war, method="GCV.Cp")
summary(gam.model9)
gam.model10 <- gam(dum_vtype2 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=binomial(link="logit"), data=boshin_war, method="GCV.Cp")
summary(gam.model10)
gam.model11 <- gam(dum_vtype3 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + s(longitude, latitude), family=binomial(link="logit"), data=boshin_war, method="GCV.Cp")
summary(gam.model11)
gam.model12 <- gam(dum_vtype4 ~ elv_avr_1000m + elv_sd_100m + Edo_100km + port_100km + highway_100km + trade_near_100km + hyde_popd_1860_100 + glues_rice_suitability_100 + avgPrecip_1901to1910_100 + avgTemp_1901to1910_100 + castle + frontline + s(longitude, latitude), family=binomial(link="logit"), data=boshin_war, method="GCV.Cp")
summary(gam.model12)

#################################################
# End of script
#################################################